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Abstract 

Two different methods are proposed for the generation of wide classes of 
exact solutions to the stationary Gross-Pitaevskii equation (GPE). The first 
method, suggested by the work by Kondrat'ev and Miller (1966), applies to one- 
dimensional (ID) GPE. It is based on the similarity between the GPE and the 
integrable Gardner equation, all solutions of the latter equation (both stationary 
and nonstationary ones) generating exact solutions to the GPE, with the poten- 
tial function proportional to the corresponding solutions. The second method 
is based on the "inverse problem" for the GPE, i.e. construction of a potential 
function which provides a desirable solution to the equation. Systematic results 
are presented for ID and 2D cases. Both methods are illustrated by a variety of 
localized solutions, including solitary vortices, for both attractive and repulsive 
nonlinearity in the GPE. The stability of the ID solutions is tested by direct 
simulations of the time-dependent GPE. 
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1 Introduction 



The Gross-Pitaevskii equation (GPE) provides for an exceptionally accurate descrip- 
tion of the dynamics of Bose-Einstein condensate (BEC) [1]. In the general case, the 
GPE is far from integrability, which was an incentive for the development of various 
methods for simulations of this equation, including finite-difference [2], split-step [3], 
and Crank-Nicolson [I] algorithms. An efficient technique for finding stationary solu- 
tions to the GPE is based on simulations of the evolution in the imaginary time [5]. A 
review of numerical methods for the GPE can be found in Ref. [6]. 

Aside from the numerical solutions, the understanding of results produced by the 
GPE requires the knowledge of its solutions in an analytical form - approximate or, if 
possible, exact. A powerful analytical method is provided by the variational approx- 
imation [7J. Another approach which simplifies the consideration reduces the three- 
dimensional (3D) GPE to an effective ID or 2D form, if the condensate is loaded, 
respectively, into a cigar-shaped or pancake-shaped trapping potential [8]. If the con- 
densate is trapped in a deep optical-lattice (OL) potential, the continual GPE may 
be further reduced to its discrete version [S] . In the case of the repulsive nonlinearity, 
the Thomas-Fermi approximation is known to be very useful [U QJJ]. The coupled- 
mode approximation is adequate for the description of settings based on double- and 
multi-well potentials [UJ. In the case when the GPE contains a rapidly oscillating time 
dependence, one may apply the averaging approximation [12] . If terms which make the 
ID GPE different from the exactly integrable nonlinear Schrodinger equation (NLSE) 
are small, one may resort to perturbation theory based on the inverse-scattering trans- 
form [T3J [14j. A number of other approximations have been developed in the context 
of the GPE, as reviewed in Ref. |15j . 

Although exact solutions of the GPE are rare, they are useful in those cases when 
they are available (see. e.g., Ref. [IS]). An example is a family of exact stationary 
periodic solutions to the ID GPE with a specially devised periodic potential, written in 
terms of elliptic functions. This analysis was performed for both cases of the repulsive 
[TT] and attractive [18] nonlinearity in the GPE, see also Ref. [19] . Exact solutions were 
also found for dark-soliton trains representing the supersonic flow of the condensate 
|2Uj . Exact localized solutions are known in the case when the nonlinearity coefficient 
is represented by a delta-function of the spatial coordinate [21] and by a symmetric 
pair of delta-functions [22]. In fact, the latter configuration provides for an exact 
solution to the spontaneous-symmetry-breaking problem. Upon a proper change of 
the notation, the GPE may be interpreted as the NLSE for spatial optical beams in 
nonlinear waveguides [23J. Accordingly, the exact solutions found in terms of the GPE 
may also find applications to nonlinear optics. 

The purpose of this work is to propose two methods for generating exact solutions 
to the GPE, together with potentials which support them. The first method is based on 
the idea proposed by Kondrat'ev and Miller [30] more than 40 years ago, namely, using 
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known solutions of nonlinear equations as potentials for other equations (a somewhat 
similar method was later proposed for the analysis of self-trapped states in nonlinear 
optics [31]). We apply it by noting that the ID time- independent GPE with the 
potential term is equivalent to a stationary Gardner equation (GE), alias an extended 
Korteweg-de Vries (KdV) equation, containing both quadratic and cubic nonlinearities, 
if the solution is proportional to the potential. Thus, one can obtain a solution to the 
GPE, along with the necessary potential, from any solution of the GE. 

The second method is based on the consideration of an inverse problem, aiming 
to construct an appropriate potential for a given wave-function ansatz representing 
an appropriate solution. The inverse problem is relevant because it may be relatively 
easy to engineer the needed trapping potential, using external magnetic and optical 
fields [24] . Recent experiments have demonstrated that, using a rapidly moving laser 
beam focused on the condensate, one can "paint" practically any desired time-average 
potential profile in ID and 2D settings [25]. A similar approach was developed in a 
different physical setting, with the objective to find models of nonlinear dynamical 
chains admitting exact solutions for traveling discrete pulses in an analytical form 
[2"o] . By means of this approach, we produce a number of novel ID and 2D stationary 
solutions. The stability of the ID solutions obtained by both above-mentioned methods 
is tested in direct simulations. 

The paper is organized as follows. In Section [21 we elaborate the Kondrat'ev-Miller 
method in the ID case and perform the stability test of the localized modes. We show 
that stationary solutions of the GPE may be constructed using not only stationary, but 
also nonstationary solutions of the GE. In the latter case, the use of the formal temporal 
variable in the GE makes it possible to obtain a wide family of stationary solutions and 
supporting potentials which depend on a continuous parameter, the relation between 
them being different from the proportionality. In the same section, we present the 
inverse method in the ID case. The stability of these solutions is tested via direct 
simulations of the time-dependent GPE. In Section EJ we construct exact solutions to 
the GPE in the 2D case (the test of their stability will be reported elsewhere). In 
particular, we construct anisotropic solutions, using the so-called lump solitons of the 
Kadomtsev-Petviashvili equation (KP1) as the respective ansatz, axisymmetric states 
with the Gaussian radial profile, and vortices with topological charges 1 and 2. The 
paper is concluded by Section 4. 

2 Exact solutions of the one-dimensional Gross— 
Pitaevskii equation 

The scaled form of the ID GPE for complex wave function (p(£,t) is well known [1]: 
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i<p t + li<p = -<p^ + u(£)<p + a\<p\ 2 <p, (1) 

where u(£) is the trapping potential, a = +1 or —1 corresponding to the repulsive 
and attractive interactions between atoms. The constant /i in Eq. ([I]) is the chemical 
potential, the objective being to find stationary solutions corresponding to a given 
value of \i. If Eq. ([I]) is derived from the underlying GPE for the 3D condensate, the 
temporal variable t and spatial coordinate £, together with function tp and normalized 
potential are related to their counterparts measured in physical units, T, X and 
as follows: t = T (ir 2 H/md 2 ), f = 7rX/d, and 



* ( * * T) = \/2^ «) CXP ("^ T " ^) • < 2 » 



1 /tt/i x 2 



where m is the atomic mass, d is a longitudinal scale determined by the axial potential, 
a s is the s-wave scattering length, while u± and R are the transverse trapping frequency 
and radial coordinate. If m is taken for 87 Rb, and d = 1.5 /im, then t = 1 and £ = 1 
correspond, in physical units, to T w 0.3 ms and X w 0.5 /mi, respectively. Finally, 
the number of atoms in the condensate is 
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2tt | | |^(X,i?,T)| 2 dX = ^^ j |^(£)| 2 d£, (4) 

— oo — oo 

where the transverse-trapping size is a± = ^h/muj^. 

2.1 Construction of stationary solutions by the Kondrat'iev- 
Miller method 

Looking for real stationary solutions to Eq. ([T]), we arrive at equation 

if" + /Mf = u(£)ip + cnp 3 , (5) 

with the prime standing for d/d£. Particular exact solutions to Eq. (|5]) can be obtained 
by way of the approach developed in Ref. [30]. To this end, parallel to Eq. (jSJ), one 
should consider the stationary GE (see [251 133] and references therein): 

< f ) "-V<f>= -a<j) 2 + a<j) 3 1 (6) 

where V and a are constants. It is well known that exact solutions to Eq. ([6]) can be 
found in terms of the Jacobi's elliptic functions. We chose one of such solutions and 
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denote it $(£)• Then, the quadratic term in Eq. fl^]) is formally factorized, and the 
first multiplier is replaced by $(£), i.e. <fi 2 = <p ■ <ft — > $(£)</>. By substituting this into 
Eq. (jnj), one obtains: 

0" - Vcf) = + o<\) 3 . (7) 

Obviously, this equation is immediately satisfied with </> = $(£). From here it follows 
that there is the one-to-one correspondence between Eqs. (J5J) and (JJJ): 

V? « — ► 0, /i ^ — ► -V, « — )• -«$(£)• (8) 

Thus, Eq. (E]) gives rise to the exact solution, <p(£) = $(£), with chemical potential 
—V, for the external potential 

u(0 = -«$(£). (9) 



2.2 Stationary solutions in the case of the repulsive nonlin- 
ear ity 

An illustrative example - the "fat soliton" of the Gardner equation. In the case of the 
repulsive atomic interactions, which corresponds to a = +1 in Eq. (0]), the approach 
outlined above may be illustrated using a particular solution to Eq. fl^]) in the form of 
so-called "fat soliton" [33]): 

$(0 = ^ [tanh (£/A + 0) - tanh (£/A - 0)] , (10) 

(9 = (1/4) In [(1 + v)/{\ - v)\ , A = 3V2/(H, V = (2/9)(av) 2 , (11) 

with free parameters a and u, the latter one taking values < v < 1. The front and 
rear slope of the fat soliton, A, depends monotonously on u, decreasing from infinity to 
A m ; n = 3^/2/ a when v varies from to 1. The width of the soliton, L, i.e. the distance 
between its front and rear segments at the half-minimum level, $(£) = Q max /2, is 




At v — > 0, the fat soliton reduces to the bell-shaped KdV soliton, whose width is 
given by L Kd y = 3\/2\n (3 + y^)/{av). In the other limit, v — > 1, it reduces to the 
"table-top soliton" (a IT-shaped mode) with L « 2A# « -3 v / 21n[(l - v)/v\/(2av). 
The minimum width, L min ps 10.1/a, is attained at v m 0.892. The local density corre- 
sponding to normalized solution fllOp . |$/a| 2 , along with the corresponding normalized 
potential, u/a 2 , in the stationary GPE for which one has $(£) as the exact solution, 
are shown in Fig. [I] for several values of free parameter v. 
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Figure 1: Normalized density |$(£)/a| 2 corresponding to solution (|T0i) (solid lines), 
and the corresponding normalized potential u{^)/a 2 (broken lines), as functions of a£. 
Lines 1 and 1' pertain to v — 0.9; lines 2 and 2' — to v — 0.999; lines 3 and 3' - to 
v = 0.99999. The horizonal dashed lines show limit values of the solution and potential 
function for the "table-top soliton" . 

The dimensionless norm of exact solution (TlOl) . which is proportional to the number 
of particles in the underlying BEC, according to Eq. (j4j), is 



As follows from here, Aid is proportional to a, while v may be treated as a free 
parameter which determines the shape of the potential and the corresponding solution. 
The chemical potential of the solution, /i = —V, is also determined by constants a and 
u, see Eq. (fTTD. 

Because wave function ffTOl) has no zeros, it may represent the ground state in the 
corresponding potential, with chemical potential /i = — (2/9) (au) 2 [see Eq. (ITT]) ]; 
whether there exist higher-order bound states with larger discrete eigenvalues of fi 
within this nonlinear problem, remains an open question. Although it is plausible 
that this solution is stable, it is relevant to test its dynamics, under the action of 
perturbations, in direct simulations of the time-dependent equation ([T]). This was 
done by means of the Yunakovsky's method [29] in a sufficiently large domain with 
periodic boundary conditions (description of the method is presented in the Appendix). 
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Examples are shown in Fig. [2] for cases when the amplitude was initially reduced or 
increased by 10% against the stationary value. 




50 100 150 200 250 

Time 



Figure 2: The time dependence of perturbed solution (fTUl) in the model with the re- 
pulsive nonlinearity, when the initial amplitude is 10% smaller than needed for the 
stationary solution. Line 1 - maxg|Re {y?(£, t)} |, line 2 - max^jlm {</?(£> t)} I- Horizon- 
tal line 3 designates the constant amplitude of stationary soliton fllOp with v = 0.99999. 
Lines 4 and 5 show the time dependence of ma.x^\ip(C,, t)\ in the cases when the soli- 
ton's amplitude was initially reduced (line 4) or increased (line 5) by 10% against the 
stationary value. 

As seen from the figure, the amplitude of the so perturbed solution varied in time 
within the same 10%, while its spatial shape was preserved. It is also seen that the 
perturbation induced oscillations between the real and imaginary parts of the solution, 
i.e. shifted its chemical potential. In the course of the simulations, the norm of the 
solution was preserved with relative accuracy ~ 10~ 7 . The latter fact attests to the 
robustness of the ground state: under the action of this sufficiently strong perturbation, 
it features no loss through emission of radiation. 

Using the above physical estimates for BEC, it is easy to estimate physical pa- 
rameters of BEC states corresponding to the fat-soliton solutions. For instance, the 
outermost configuration in Fig. [1] represents to the potential well of the depth ~ 1 
recoil energy corresponding to d — 1.5 fim, and width L ~ 30, which is ~ 15 /xm in 
physical units (see above). These values are quite realistic for the experiment [Tl l2"4l I2"5] . 
Further, taking an experimentally relevant values of a s = 5 nm and a_i_ = 3 /im, Eq. (J4]) 
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yields the largest number of 8 Rb atoms which may form the fat soliton, N ~ 30, 000. 
This estimate shows that the constructed soliton solution is quite relevant to the ex- 
periment. Finally, the period of large-amplitude oscillations of the perturbed stable 
state, shown in Fig. [2], is ~ 100 ms. 

Other stationary solutions related to the Gardner equation. To consider various 
exact solutions to Eq. one may write the equation in the "energy conservation" 
form, (X') 2 /2 + P(X) = E, where X = (f)/a, E is a constant of integration, and P(X) 
is the effective potential, 

P+(X) = - (a/4) a 2 X 2 [(X - 2a/3) 2 + 2aW - 4/9] , (13) 

with W = V/a 2 . It is shown in Fig. [3] for a = 1. 

kP(X)/a> 




Figure 3: Normalized polynomial P(X)/a 2 defined in Eq. f[T3"|) with a — 1 at different 
values of W: line 1 - for W = 9/32, line 2 - for W = 1/4, line 3 - for W = 17/72, line 
4 - for W = 2/9, line 5 - for W = 1 /5, line 6 - for W = 0, and line 7 - for W = -2/9. 

For W > 1/4, the polynomial has a single real maximum at X = 0, hence neither 
periodic nor solitary solutions are possible. At W < 1/4, it has three real extrema, at 
points X = 0, (1/2)(1 - y/l-AW), (1/2)(1 + Vl -4W). In that case > first appears 
a depression-type solitary solution, in the form of a "bubble" against the constant 

background value of the field, 0o = ( a /2) ^1 + a/1 — AV/a 2 ^ . The typical potential 
profile corresponding to the "bubble" is depicted by the line 3 in Fig. [3j 

At W = 2/9, the potential becomes symmetric, as shown in Fig. E] by line 4. 
Solution ( TTUj) with u = 1 corresponds exactly to this case. For smaller values of W, 
when it varies from 2/9 to 0, the right maximum of the potential function becomes taller 
then the left one, making it possible to have bright solitons with the zero background. 
They correspond to solution ( fTOj) with v < 1. The solution vanishes at W — > +0, 
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i.e. v —> +0. For W < 0, the left maximum of the potential shifts from the origin to 
the left, see Fig. E}. In this case, solitons exist against the negative background, with 

0o= («/2) (l-y/l-W/a*y 

Periodic solutions of Eq. (JBJ) can be analyzed similarly. The corresponding station- 
ary solutions of the GPE can be readily obtained by means of the method described 
above, in terms of elliptic functions, similar to the periodic solutions reported in Ref. 



2.3 Stationary solutions in the case of the attractive nonlin- 
ear ity 

In the case of the attractive nonlinearity, i.e. a = — 1 in Eq. (j5J), the potential function 
is different from that shown in Fig. [3J It is shown in Fig. [5] in two different scales, as it 
is impossible to display all details using a single scale. For W > —1/4, three possible 
real extrema of this polynomial are located at points 



x = o, ~ (l - v / TT4f) , ~ (l + VT 



AW 



otherwise the polynomial has a single real minimum at X = 0, hence solitary solutions 
do not exist for W < —1/4. 



kP(X)/a 2 




Ap(X)/a 2 




a) b) 

Figure 4: Normalized polynomial P_(X)/a 2 defined by Eq. (I13|) with a = — 1 for 
different values of W. Frame (a): line 1 - for W = —0.26, line 2 - for W = —0.25, line 
3 - for W = -0.23, line 4 - for W = -0.22, and line 5 - for W = -0.21. Frame (b) 
uses a different scale: line 6 - for W = —0.22 [the same as line 4 in (a)], line 7 - for 
W = -0.11, line 8 - for W = 0.11, and line 9 - for W = 0.44. 
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The first solitary-type solution emerges at W — — 1/4. In this case, the potential 
still has only one minimum at X = 0, but the inflexion point appears at X = —1/2 (the 
corresponding potential function is shown in Fig. H]by line 2). A particular solution 
corresponding to W = —1/4 represents the algebraic soliton sitting on top of a pedestal 
(constant- value background), as a solution to Eq. ([5]) with free parameter a and V = 
-fi = -a 2 /4: 

*«> = f {Irdm - ') ■ (14 » 

At W > —1/4, one more minimum appears in the potential profile (see, e.g., line 3 
in Fig. |4]). In this case, two families of solitons on a pedestal are generated by Eq. ([6]): 

* ± tt) = a J 1 3(1 + 2*) 1 (15) 

[ l + 3r/± v/1 + 3z//2cosh af + 2i/) J 

where v is a free parameter ranging between —1/2 and 0, and V = — \i = o?v{\ + v). 
The local densities corresponding to solution $+(0 has the form of a double dark 
soliton (with two zeros), whereas solution $-(0 is shaped as a bump on top of the 
pedestal. The local densities corresponding to solutions (fl~5]) . along with the respective 
potentials, for which these are exact solutions to the GPE, are shown in Fig. [5] for 
v = -5/12. 

When W increases further and becomes equal to —2/9, potential function f)13p 
becomes symmetric with respect to the vertical line X = —1/3 (see line 4 in Fig. HJ), 
getting then asymmetric, with the left minimum falling deeper than the right one when 
W increases further (see, e.g., line 5 in Fig. H]). For the particular case of W = —2/9, 
solutions ffT5"j) reduce to 



1 ± v^sechf 

V\/3 



For W ranging between —2/9 and —1/9, the potential function ( TT3|) becomes again 
asymmetric, as mentioned above, while the corresponding solutions are still given by 
Eq. f)15p with —1/3 < v < 0. They remain unstable, because of the nonzero back- 
ground. 

At W — — 1/9, the maximum and minimum of the potential merge at X = 0. An- 
other inflexion point emerges in the potential profile in this case [see line 7 in Fig. IHb)]. 
The corresponding solution to Eq. §6§ is an algebraic soliton with V = —fi = and 
zero background, 

{q) 3 l + 2a 2 e/9 

However, as follows from Eq. ([9]), this solution corresponds to the maximum of the 
physical potential, hence its is unstable (which was confirmed by direct simulations). 
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Figure 5: Solutions (|T5|) in terms of ($±/a) 2 (solid lines), and respective normalized 
potentials u/a 2 (broken lines) are shown as functions of normalized coordinate a£. 
Lines 1, 1' and 2, 2' pertain, respectively, to signs plus and minus in Eq. (j!5p . In both 
cases, v = —5/12. 



For W > —1/9, two families of exponentially localized solitons are generated by Eq. 

*±(0 = — ; > 16 

1± v/l + 9z/ 2 /2cosh(^0 



where free parameter v > determines the inverse width of the soliton, and V = —fi = 
{civ) 2 . These solutions and corresponding potentials are shown in Fig. [6] for v = 1. As 
follows from Eq. solution is unstable, as it represents a soliton sitting at 

the potential maximum. However, solution $+(0 is trapped in the minimum of the 
attractive potential, hence it may be stable. At v — > 0, the latter solution smoothly 
vanishes, whereas the unstable one reduces to the above-mentioned unstable algebraic 
soliton. 

The stability of all the solutions found in the model with the attractive nonlinearity 
was tested via direct simulations of Eq. ([T]). First, the expected instability of the 
solutions on the pedestal and localized modes placed at the maximum of the potential 
was corroborated. In the former case, the modulational instability of the background 
leads to the formation of a chaotic "gas" of interacting solitons. An examples of 
such instability onset is shown in Fig. [71 which shows a spatial distribution of the 
density, | | 2 , at different moments of time when the spatial period is L = 128. 
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Figure 6: ($±/a) 2 for solutions (fl6|) with v — 1 (solid lines), and the corresponding 
normalized potentials, u/ot 2 (broken lines), as functions of a£. Lines 1, 1' and 2, 2' 
pertain, respectively, to signs + and — in Eq. (1T5]1 . 

In the latter case, small random perturbations may either cause the soliton to roll 
down from the unstable position (see an example in Fig. [Sk), or split - symmetrically 
(see Fig. [Hb), or sometimes asymmetrically (see Fig. [5b) - into two solitons moving 
in opposite directions. In particular, the splitting was, naturally, observed under the 
action of an initial perturbation which made the amplitude of the unstably pinned 
soliton smaller, hence making it more similar to a quasi-linear wave packet which is 
subject to the splitting by the potential barrier. In fact, the strongly asymmetric 
splitting may be realized as a result of strong emission of radiation from the unstable 
soliton and self-retrapping of the emitted wave packet into a small-amplitude soliton. 
Because the simulations were run in the domain with periodic boundary conditions, 
we also observed that, in the case of the splitting of the unstable soliton into two, like 
in Figs. [Hb and [Bb, the secondary solitons survived in the course of numerous head-on 
collisions in the course of their circular motion. If the instability gave rise to a moving 
soliton and a radiation wave-train, the secondary interactions between them would not 
destroy the soliton either. 

Also in agreement with the expectation formulated above, the numerical solutions 
demonstrate the stability of solitons $+ gives by Eq. (TT6]) . In this case, the amplitude of 
the perturbed soliton, max^|<^(^)|, varies in time around some average value, as shown 
by lines 4 and 5 in Fig. [9], and remains within the same range of the deviation from the 
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Figure 7: An example of the instability development for the algebraic soliton on a 
pedestal, given by Eq. ( f!4|) in the model with the self- attraction: (a) - t = 0; (b) - 
t = 30; (c) — t — 40; (d) - t = 80; (e) - t = 100 (the vertical scale in the three latter 
cases is compressed by a factor of 2.5). 
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Figure 8: The evolution of a disturbed exponentially localized soliton from Eq. 
(fT5"j) (only a part of the total spatial period of L = 128 is shown). Panel (a) - the 
initial amplitude is 10% greater than the stationary value corresponding to a — v — 1 
(line 1 - t = 0, line 2 - t = 12, line 3 - £ = 16, line 4 - * = 18, line 5 - t = 20, and 
line 6 - £ = 22). Panel (b) - initial amplitude is 10% smaller than the stationary value 
(line 1 - t = 0, line 2 - i = 2, line 3 - t = 6, line 4 - i = 10, and line 5 - 1 = 14). Panel 
(c) shows the evolution of initially undisturbed soliton (fT6|) . with a = 1 and v = 0.675, 
under the influence of small errors of the numerical truncation (line 1 - t = 0, line 2 - 
t = 8, line 3 - t = 10, line 4 - i = 12, and line 5 - t = 14). 
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stationary value as the initial perturbation. The stability of these solitons is similar 
to that which was demonstrated above for the table-top soliton in Fig. [5] in the model 
with repulsive nonlinearity. 




50 100 150 200 250 
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Figure 9: The evolution of perturbed solution $ + from Eq. (fl6]) in the model with 
the attractive nonlinearity, when the amplitude of the initial perturbation was 10% 
smaller than needed for the stationary solution. Line 1 - max^l Re t)} |; line 2 
- max^llm {</9 (£,£)} | . Horizontal line 3 shows the constant amplitude of stationary 
soliton $ + with v — 1; lines 4 and 5 show the time dependence of max^l £)| in cases 
when the soliton's amplitude was reduced (line 4) or increased (line 5) by 10% against 
the stationary value. 

2.4 Other stationary solutions obtained from nonstationary 
solutions of the auxiliary Gardner equation 

The nonstationary version of GE (jSJ), 

<P T + (c0 + -<j?- a0 3 ) € = 0. (17) 

can also be used for the purpose of generating stable stationary solutions to the GPE 
[note that here, in comparison with Eq. (JBJ), we set a = — 1]. We stress that formal 
temporal variable r in this equation has nothing to do with physical time t in Eq. (1T1). 
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Equation (|17|) is tantamount to the integrable modified KdV equation 
, It may be formally integrated once in £ and cast in the following form: 



c6 + 



d 2 



crc 







$(£,*) dr 



<b(x', t) dx' 



where <&(x,t) is one of the particular nonstationary solutions to Eq. (I17j) . The term 
on the right-hand side of Eq. (fT8j) can be combined with term $(£, r)0 on the left-hand 
side, to generate the stationary GPE in the form of Eq. (jSj) with c — \i and potential 
' u (£) r ) which depends on free parameter r: 



u(£,r) = $(£,r) 
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<b(x', t) dx', 



(19) 



cf. Eq. for the case when the stationary GE was used. It is worthy to stress 
that, in the present case, the solution and the potential which support it are no longer 
proportional to each other. 

As an example, we take, following Ref. [35], a solution to nonstationary GE f|T7|) 
with a = 1, which describes disintegration of the initial configuration into two fat 
solitons: 



$(£,r) 



v\ - v\( 1 1 



where 



Z±± = V\ tanh 



Z2± = V2 coth 
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(20) 



(21) 



(22) 



81,2 = i^)' 1 In [(1 + 1/1,2) / (1 - ^1,2)] • (23) 

Two examples of these solutions are shown in Fig. [10] for a fixed value of V\ = 0.75 
and two different values of the other parameter, v-i = 0.8 and 0.999. The configuration 
shown in Fig. IToTb) actually represents a pair of table-top Gardner solitons ffTUj) with 
different parameters. Solution (I2"01 - fl2lll) may be treated as a continuous family (pa- 
rameterized by r) of stationary solutions to Eq. (jSJ), with the corresponding potential 
produced by Eqs. (I20l-(l23l . 

To conclude this subsection, it is relevant to mention exact nonstationary solutions 
in the form of breathers, for Gardner equation (ITTj) corresponding to the GPE with 
the attractive nonlinearity (a = —1) [34] . Such breathers look as two periodically 
interacting exponentially localized solitons, or as envelope solitons of the NLSE. At 
any fixed value of r in Eq. ()17[) . the breathers generate, as outlined above, exact 
stationary solutions of the GPE with the corresponding potentials given by Eq. (Il9j) . 
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Figure 10: Solutions to the nonstationary Gardner equation (II 7p with c = and (7 = 1, 
which describe disintegration of the initial configuration into a pair of "fat solitons" . 
They emerge with parameters ui = 0.75, v<i = 0.8 in (a), and v\ = 0.75, z/2 = 0.999 in 
(b). Lines 1, 2, and 3 correspond, respectively, to r = 0, r = 1500, and r = 3000. 

2.5 Reconstruction of the supporting potential in the GPE 
for an arbitrary matter-wave distribution 

An arbitrary distribution of the stationary matter- wave field, y?(£); can be made an 
exact solution to stationary GPE ([5]), if the potential in the equation is chosen as 

«(0 = " V(0 + A*- (24) 

Below, we demonstrate that this, seemingly "trivial", approach may also produce es- 
sential results. 

An example: the Gaussian profile of the matter wave. First, we take a Gaus- 
sian matter-wave pulse, which is the case of obvious interest to applications, = 
Aexp [— (£//) 2 ]. Substituting this into Eq. f )2i|) . one finds: 

«(£) = I (2^ " l) " °^ exp [-2(i/lf] + im. 
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The trial solution, (p(£), and the corresponding potential can be presented in the di- 
mensionless form: 

F(C) = exp(-C 2 ), 1^0 = Jj(2C 2 -l)-^exp(-2C 2 ) + M, (25) 

where C = F{() = <p(()/A, v e (() = u{Q/A\ S = Al, M = fi/A 2 . While F(() does 
not contain any parameter, the dimensionless potential depends on two independent 
constants, S and M, for each sign of a = ±1. Figure [TT] shows the squared normalized 
solution, F 2 ((), and the corresponding potentials for both signs of a, as given by 
Eq. fl25|) for several values of S and M — 0. Note that, in the case displayed in 
Fig. [TTb . the potential corresponding to the GPE with the self-attraction nonlinearity 
(a = —1) features a double-well shape. As follows from Eq. 025p . such a shape may 
occurs only in the case of a — — 1, provided that S exceeds a threshold value, S t h T = y/2. 

A derivative-Gaussian profile of the matter wave. Trial function </?(£) used above is 
an even function of £ without nodes, which, apparently, represents the ground state for 
the nonlinear GPE with the given potential. Here we aim to consider another example, 
when the trial function is chosen as an odd one, with a single node, thus representing 
the first excited state. To this end, we take = A£ exp [— (£//) 2 ] and derive the 
corresponding potential from Eq. ( |24|) : 

«(£) = | (2^ - 3) - aA 2 e exp [-2(i/lf] + im. 

It is again convenient to present the trial solution, <p(£), and the corresponding potential 
in the dimensionless form: 

F(C) = -Cexp (-C 2 ), v e (C) = A (2C 2 - 3) - aC 2 exp (-2( 2 ) + M, (26) 

where, this time, C = 0, F(Q = <p(Q/A, u e (C) = u(()/A 2 l 4 , S = Al\ and M = 
fj,/(Al) 2 . Plots corresponding to this trial function and the supporting potentials are 
displayed in Fig. [T2]for M = and both signs of the parameter a. Apparently, function 
F(() represents the first excited eigenmode in the corresponding potential well v e ((). 
Note that, as it follows from Eq. (126|) . the trapping potential has a single-well structure 
at S < S^. = 2 for a = +1, and at S < S^ T = 2e for a = —1. In the opposite case, 
the potential acquires the double- well structure when S > S^ r for a = +1 (see line 
1' in Fig. IT2"b). and the triple- well structure when S > S^ r for a = — 1 (see line 2' in 
Fig. [T2"b). The shapes of the potential at the critical values of S = S^ v are shown in 
Fig. d2b. 

A comb-top-Gaussian profile of the mater wave. Here we consider the trial function 
in the form of the Gaussian with a superimposed "comb" , which corresponds to the 
physically relevant combination of an OL and external parabolic trap: 

<p(£) = A exp {-e/l\) + B cos (kx) exp (-£ 2 // 2 ), (27) 
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-0.1 



FIG. 11 FIG. 12 

Figure 11: Normalized Gaussian solution f l23|) in terms of F 2 (() (solid curve 1) and 
corresponding normalized potentials f e (C) (broken curves 1' and 2', pertaining to a = 1 
and a — — 1, respectively), as functions of normalized coordinate (. Panels a), b), and 
c) were generated for S — 0.5, S = 1, and S = 5, respectively. Note that the plots are 
shown on different scales. 

Figure 12: The same as in Fig. [11] but for the derivative-Gaussian solution (126]) . The 
corresponding potentials are additionally reduced by the factor K: v e (()/K, where K 
is different in each panel. Panel a): S = 1, K — 100 for a = ±1; panel b): S = 2, 
K = 20 for a = 1 and S = 2e, K = 4 for a = -1; and panel c): S = 25, K = 2 for 
cr = ±1. 
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where A, B and l\^ are arbitrary constants. This function resembles, in particular, a 
numerical solution which was found in Ref. [H] (see also review [5]) for the OL potential. 
Function (1271) and the corresponding potential, as given by Eq. ( 1241) . can be presented 
in the following dimensionless form: 

F(C) = exp (-C 2 ) + 6 cos (O exp [-(s() 2 }, (28) 



Ve(C) 



£ 4 C 2 - T " T cos («0 + £ < sin «) 



exp [— [e 



S 2 l + 6cos(KC)exp[-(e 2 - 1)C 5 

- a exp (-2C 2 ) {1 + 6 cos «) exp [-(e 2 - 1)( 2 ] } 2 - 



(29) 



where F{Q = <p({)/A, v e (£) = u(i)/A\ and C = £/h, b = B/A, k = kl x , e = h/l 2 , 
S = Ah, M = —fi/A 2 . Varying parameters S, M, k, b and e, one can obtain a wide 
class of solutions. The corresponding potentials asymptotically approach the parabolic 
shape at large |£|, featuring a complex oscillatory shape at the center. Solution ( 128]) . 
F 2 (() and the corresponding potential (1291) are shown in Fig. [T3]for both signs of a. 



The stability of the Gaussian-type solutions. Stability of all solutions presented in 
this section was tested via simulations of Eq. ([1]). The results are summarized as 
follows. 

(1) In the case of the repulsive nonlinearity, a — 1, the Gaussian solution with 
potential (1231) is stable for all values of S. It is stable too in the case of a — —1 (the 
attractive nonlinearity) if the corresponding potential features the single-well shape 
(such as shown in Figs. ITTk and 9b), i.e. S < SW = \/2 5 see above. However, in the 
model with the attractive nonlinearity, the solution naturally becomes unstable when 
the single-well potential transforms into the double-well potential, i.e. when S > y/2 
(see line 2' in Fig. ITTb). In the latter case, the solution preserves its shape until t < 20, 
and then spontaneously splits into two pulses which quasi-regularly oscillate relative 
to each other. 

(2) The derivative- Gaussian solution with potential ( )26|) is stable in both cases of 
cr = ±1, provided that the underlying potential keeps the single-well shape, i.e. until 
S < S+ r = 2 for a = +1 and S < = 2e for a = -1 (see Figs. HSi and DJd). At 
greater values of S, the solutions in the double-well potential (see Fig.[T2"b) are unstable, 
for either sign of a . However, manifestations of the instability are different for o = +1 
and a = — 1. In the former case, the initial distribution was preserved in a quasi-stable 
state until t 40. Then, the profile of became asymmetric, with one maximum 
(for instance, the right one, as in Fig. [T4"k ) being greater than the other. After reaching 
a well-pronounced asymmetric shape, the process reverted, making the left maximum 
greater than the right one. This process repeated persistently, so that the initial soliton 
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2.0 A 




Figure 13: Squared solution f[2"51) . F 2 ((), (panel a) and the corresponding potential fl2"§|) 
(reduced by the factor of 10) for a = 1 (panel b) and a — — 1 (panel c) as functions of 
C Parameters are S = 5, M = 0, k = 40, 6 = 0.25, e = 2. 
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Figure 14: The spatial profile of \ip\ 2 at different times, for the derivative- Gaussian 
initial configuration with S = 25 (the size of the computation domain is L = 128). In 
panel (a) pertaining to the case of a — 1 lines 1, 2 and 3 correspond, respectively, to 
t = 0, 50 and 60. In panel (b) pertaining to the case of a = — 1 lines 1, 2, 3, 4 and 
5 correspond, respectively, to t = 0, 1, 2, 3, and 4 (note that there are two pulses of 
smaller amplitudes near the center at t = 1, 3 and 4). 

was eventually transformed into an immobile breather consisting of two non- stationary 
pulses which oscillate in time quasi-randomly due to the energy exchange between them 
and, apparently, due to their interaction with a linear wave-train shed off by the pulses. 
It is worthy to note that the instability of the odd mode trapped in the double-well 
potential shown, for instance, in Fig. fT2b. sets in via the breaking of the skew symmetry 
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of this mode, in agreement with the general principle that the repulsive nonlinearity 
gives rise to the symmetry breaking of odd modes trapped in double-well potentials 



In the case of a — — 1, the character of instability observed at S > is different. 
For instance, at S = 25 the solution became unstable at £ ~ 5, splitting into four pulses 
and a small-amplitude wave train (see line 2 in Fig. lT4"|) . Two of those four pulses moved 
to the left and two others - to the right. After passing some distance, these pulses 
bounced back from the parabolic potential, moved back towards the center, and formed 
two very narrow and closely located spikes (see line 3 in Fig. HI]) . Such cycles of the 
splitting and partial recombination with the pulse compression in the vicinity of the 
center repeated indefinitely long. Thus, some type of a double breather is formed in 
this case too. 

(3) The comb-top Gaussian solution with potential (T2l?j) is stable in the model with 
the repulsive nonlinearity, a = 1 (we did not explore the stability of this solution in 
a wide range of parameters; here we report an example for 6 = 1/4, e = 2, k = 40, 
M — 0, and S = 5). Namely, if external perturbations were added to the solution, 
e.g., A in Eq. (J27J) was taken 10% smaller or grater against its stationary value, the 
evolution lead to time variations of \ip\ within the same 10% range, similar to what was 
reported in item (1) above for the Gaussian initial distribution. 

In the case of the attractive nonlinearity, a — — 1, with the same set of parameters 
as above, the initial real wave function (12 7|) was quickly, within t « 5, transformed 
into a complex one, with a subsequent quasi-random energy exchange between the 
real and imaginary parts. At the initial stage of the evolution, the central part of the 
solution would transform into a very narrow large-amplitude pulse, with two side wings 
represented by small-amplitude wide pulses. Then, the central pulse would oscillate in 
time quasi-randomly. 

3 Exact solutions for the two-dimensional station- 
ary Gross— Pitaevskii equation 

In this section we proceed to the consideration of the 2D version of the stationary GPE, 
taken in the dimensionless form similar to Eq. (j5J), with spatial coordinates £ and r\: 



The approach similar to that which was developed in subsection I2.5l can be employed to 
construct an appropriate 2D potential on the basis of a given solution. From Eq. (130]) 
one formally deduces 



nu. 



<P& + <Prm = [«(£ , rj)-tA<p + a M V 



(30) 



(31) 
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Examples presented below aim to demonstrate solutions which may be useful to phys- 
ical applications. 

The Kadomtsev-Petviashvili lump soliton. As the first trial function, we take an 
anisotropic 2D weakly localized ansatz known as the lump solution to the KP1 equation 

m 3 - {j/af + ( V /b) 2 

[3 + (e/«) 2 + (^) s 

Substituting it as <p into Eq. (13T|) . one obtains the corresponding potential, 



Htn) = i*A ,:-)y:>^ )■",:> (32) 



6 PAP rfi 3 - ^ + ?] 

= -^T^^SS^T - 144a 7 2 y ~ + M ' 

(3 + e +^ 2 ) 2 (3-e + ?f) (s+f+v 2 



2 

■2 



where w e = u/A 2 , £ = £/a, rj = rj/b, S = Ab, (3 = b/a, M = p/A 2 , and P^rj) = 
tf* - l)(f + rf) + 6(1 - /3 2 )(^) 2 - 2(f + rj 2 ) ~ 9/? 2 (2? + 1) + 3(1 + 3(3 2 ). In the 
particular case of (3 = 1 (a = 6), Eq. (TH3T) simplifies to 

= _ i4) 2 ,3- 5 r,,, r ,3 r +ff + M 

(3 + e + ^ 2 ) 4 

Both expressions f[3"3"]) and (|34p correspond to anisotropic 2D trapping potentials. Fig- 
ure ITo'a shows a 3D view of lump solution (I3"2"|) for a = b, and Figs. [ToTj.c display the 
corresponding potential (154"]) for a = ±1. In the case a = — 1, the potential represents 
a 2D hump, i.e. it is repulsive, hence the corresponding solution ( 132]) is apparently 
unstable. 

2D Gaussian trial function. Another natural example in the 2D case is provided by 
the solution ansatz in the form of an axisymmetric Gaussian, $(r) = Aexp (— r 2 // 2 ), 
where r 2 = £ 2 + rf. From Eq. f |30|) we deduce the potential in the dimensionless form: 

<P) = ^(P 2 ~ 1 )- ffex P (- 2 P 2 ) + M > ( 35 ) 

where v e = u/A 2 , p = r/l, S = Al and M = fi/A 2 (the same particular solution was 
recently obtained in Ref. [36] in a different way, and its stability has been established 
in direct simulation). The principal cross-sections of the squared Gaussian solution 
and the corresponding potentials for a — ±1 are shown in Fig. [16] for three values of 
S. In the case of the attractive nonlinearity, a = — 1, potential function v e may feature 
a local maximum at the center, which appears at S > v2 (see Fig. IT5"b). 
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Figure 15: (a) Lump solution (|32|) shown in terms of [$(£,77) /A} 2 for a = b and 
M = 0; and the corresponding potential (IMj) for er = 1 (b) and a = — 1 (c). 
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Other examples represent 2D patterns of a different type, namely, vortices. First 
we present the unitary vortex, with topological charge J = 1 (it may be treated as 
a 2D counterpart of the ID derivative-Gaussian profile considered above): $(£,77) = 
A(£ + irj) exp (— r 2 // 2 ), where r = a/£ 2 + rj 2 . For the complex solution it is convenient 
to cast the stationary GPE flU into the following real form: (1/2) [A (|<p| 2 ) - 2|V<p| 2 ]- 
o~\(p\ A = 77) — p]|<p| 2 , where A is the 2D Laplacian. From this equation, one can 
deduce the potential: 



A 


<fi 


2 


1 


2|V^| 2 


2 




2 



Substituting here the wave function of the vortex, the potential can be obtained in the 
explicit form [cf. Eq. (12 6]) ]: 

v e (p) = (4/S 2 ) (p 2 - 2) - ap 2 exp (-2p 2 ) + M, (37) 

where w e = u/(Al) 2 , p — r/l, S — Al 2 , and M = fi/(Al) 2 . Principal cross-sections 
of the solution for |$(p)| 2 and the corresponding potentials for a — ±1 are shown in 
Fig. El for three values of S. 

The double vortex. We have also considered the trial solution in the form of the 
vortex with J = 2, viz., $(£,77) = A(x + iy) 2 exp (— r 2 / 7 2 ) . One can readily deduce 
from Eq. the potential supporting this solution: 

v e (p) = (4/S 2 ) (p 2 - 3) - ap 4 exp (-2p 2 ) + M, (38) 

where v e = u/(Al 2 ) 2 , p = r/l, S = Al 3 , and M = fi/(Al 2 ) 2 . Principal cross-sections 
of the solution for |$(p)| 2 and the corresponding potentials for a = ±1 are shown in 
Fig. US] for the same three values of 5" as in Fig. [T71 



4 Exact solutions for the three-dimensional station- 
ary Gross— Pitaevskii equation 

In this section we briefly demonstrate the applicability of the above method to the 3D 
case (with the usual Laplacian corresponding to the GPE). The scaled 3D version of 
the stationary GPE is 

d 2 (p d 2 tp d 2 f v / j- \ 1 ,9 / x 

0^+0^+ d(^ = *) + d V + ff M V, (39) 

where function depends on three spatial coordinates - £, 77, and £. From Eq. flH^j) 
one can readily deduce the necessary potential: 
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a) 




c) 



Figure 16: The radial distribution of density |y?(p)| 2 for the 2D Gaussian solution (lines 
1 in each panel) and the corresponding potentials (1351) with M = for the repulsive 
(a = 1, lines 1') and attractive (er = —1, lines 2') nonlinearities: a) S = 1, b) S = v2> 
c) S = 3. 
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Figure 17: The radial distribution of density |<y9(p)| 2 in the unitary vortex solution 
(lines 1), and the corresponding potentials ( |37|) with M = for the repulsive (a = 1, 
lines 1') and attractive (a = —1, lines 2') nonlinearity. a) S = 3, b) S = 4.5, c) S = 10. 

Figure 18: The radial distribution of density |y9(p)| 2 for the double vortex (lines 1), 
and the corresponding potentials (1381) with M = for the repulsive (er = 1, lines 1') 
and attractive (er = —1, lines 2') nonlinearities. a) 5* = 3, b) S — 4.5, c) S = 10. 
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_ , 1 (d 2 <p d 2 y d 2 <p\ 2 

«™ = t{w + w + w)-' v (40) 

We here present an example of the 3D trial function in the form of the spherically 
symmetric Gaussian, $(r) = Aexp (— r 2 /l 2 ), where r 2 = £ 2 +7/ 2 +C 2 - The corresponding 
potential can be found from Eq. ( |39|) : 

<P) = ^ (p 2 " " o- exp (-2p 2 ) - M, (41) 

where v e = u/A 2 , p = r/l, S = Al, and M = p/A 2 . The plot of the potential, v e (p), is 
similar to that shown in Fig. [16j In a similar way, one can readily construct potentials 
for many other cases including 3D generalizations of the examples presented above in 
terms of the 2D geometry. 

5 Conclusions 

We have demonstrated that numerous exact ID stationary solutions to the GPE 
(Gross-Pitaevskii equation) may be constructed with the help of the Kondrat'iev- 
Miller method [3U] . Within the framework of this method, the corresponding potential 
function in the GPE is proportional to stationary solution y2(£), which, in turn, was 
taken as a solution to the stationary GE (Gardner equation). The stability of the ID 
solutions was tested through direct simulations of the time-dependent GPE. It was 
found that some solitary-type solutions are stable - in particular, those corresponding 
to the solution of the GE in the form of the "fat" (table-top) soliton, which is given by 
Eq. fllOp . in the case of the repulsive nonlinearity. A stable soliton solution in the case 
of the attractive nonlinearity was also found, viz., the exponentially localized solution 
given by expression ffT6]) . 

Further, we have proposed an "inverse method" for the GPE, as a way to construct 
appropriate potentials for a given distribution of the wave function. It was demon- 
strated that this method helps to produce many solutions in ID and 2D settings. The 
stability of all the so found ID solutions has been tested in direct simulations. The ID 
and 2D potentials constructed here as supports for basic natural types of the localized 
matter-wave distributions are fairly simple, and may be realized in the experiment by 
means of currently available techniques, based on the design of appropriate magnetic 
and optical traps for the BEG 

The numerical scheme employed here for the simulations of the time-dependent 
GPE in ID is based on the Yunakovsky's method of the operator exponential, which 
has been used in many previous works (see, e.g., Ref. [29]). The method is also efficient 
in obtaining solutions to NLSE and GPE in the space of any dimension. It is briefly 
described in Appendix. 



29 



In addition to testing the stability of the 2D localized solutions, another remaining 
issue is to check whether the specially designed potentials which support ID and 2D 
solutions, as ground states, may also sustain higher-order bound states in the respec- 
tive nonlinear Further, the inverse method can be readily extended to the 3D settings. 
Some results have been already obtained in this direction, to be reported elsewhere. 
Finally, it may be quite interesting to apply both the GE and the inverse method to 
constructing exact solutions for dark solitons in ID and circular dark solitons in 2D 
[37] . in the case of the modulationally stable background. These generalizations will 
be also reported elsewhere. 

Acknowledgements. This work was partially supported by the German-Israel 
Foundation through the grant No. 149/2006. Y.S. appreciates hospitality of the Faculty 
of Engineering at the Tel Aviv University during his visit in 2009. 

6 Appendix: The numerical algorithm for simula- 
tions of the Gross-Pitaevskii equation 

In this Appendix we describe a numerical algorithm based on the method originally 
developed by Yunakovsky for the numerical solution of the NLSE [29]. The method 
works equally well in ID, 2D and 3D settings. Below give a brief account of the 
method in the application to the ID case, it generalization for 2D and 3D settings 
being straightforward. 

One starts by the application of the Fourier transform to variable £ in Eq. (JT]): 

= k 2 $ + F{[u(0 + fi]<p} + aF {\<p\ V} , (42) 

where k is the respective wavenumber, and the tilde stands for the Fourier image, which 
is generated by the Fourier-transform operator F. Next, we introduce a new function, 
u(k,t) = <j5(k,t)exp (ik 2 t), and rewrite Eq. fflS]) accordingly: 

i— = F { [U^C) + v\ V \ 2 } if} exp (ik 2 t), (43) 

where t/ M (£) = u(£) + ji. This equation may be formally integrated in t, yielding, in 
terms of tp, 

t 

fi(k, t) = fi(k, 0) exp (-ikH) - i J F{ [U^) + a\(p\ 2 ] ^} exp [-ik 2 {t - t)\ dr. (44) 

o 
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The integral on the right-hand side can be approximately calculated, over a small time 
interval, with the help of the trapezoidal rule. The result, valid up to 0(t 3 ), is 



x 



ip(k,t) = (fo exp (— ik 2 t) — (it/2) 
F { [Up(£) + a|^o| 2 ] ^0} exp {-ikH) + F { [U^) + a\ip\ 2 ] if) 



(45) 



Next, we collect on the left-hand side those terms which depend on the current time, 
and leave on the right-hand side the terms which depend on initial conditions: 



it 



<p(k,t) + -F {[U&) + a\<p\*\ = 
<?o - % \f { [U^iC) + a\<p \ 2 ] <Po} ] exp (-ik 2 t). 



By applying the inverse Fourier transform, F 1 , to Eq. (146|) . one obtains 
it 



(46) 



[u»(0 + °\^t)\ 2 ]\^t) 



F~ l 



it 



<p (k) - -F {[11,(0 + a\ Vo \ 2 } Vo } 



exp(-ikH) } = B(£,t). (47) 



Function B(£,t) produced by Eq. (14"7|) is an explicit result obtained from the given 
initial condition, <po(k) = <p(£,0). Then, function </?(£, i) can be formally found from 
Eq. flED: 

B(£ t) 

^ t] = 1 + W/2) WM+^WWY (48) 

To make this formula practical, one needs to define \ip(C,,t)\ 2 in the denominator of 



Eq. (1451) . To do that, take Eq. ([47]) and multiply it by the complex conjugate counter- 
part, which yields 



1 + (t 2 /4) [u,(o +a\ip\ 2 ] 2 \ we, t)\ 2 = m, t) 1 



(49) 



Thus, denoting z = \(p(£,,t)\ 2 and taking into account that a 2 = 1, we obtain a cubic 
equation for z: 



+ 2aU,(i)z 2 + [(4/t 2 ) + U,(0] z - (4/t 2 ) |£(£, t)\ 2 = 0. 



(50) 



The cubic equation can be solved analytically, in principle. Its single real root (two 
others are complex) can be found by means of symbolic calculations realized by means 
of software such as Maple: 



2 

z = — 
3 



Root3 



l/t 2 + ^(0/4-^(0/3 
Root 3 



^(0 



(51) 
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where we define 



Root3 



X 



i 

T 2 



+ 



27 

it* 



1/3 



3 {[1 - £/,(£)] Ufc)t« + 4 [3 - 2^(0 + a (9 - 8t^(0) l^(e,t)| 2 ] Ufc)t* + 
[48^(0 + (144cr^(0 + 108) \B(^t)\ 2 - lQ8aU 2 (0] t 2 + 64} . (52) 



Once root z was found, it can be substituted into the denominator of Eq. (14"5|) ; then, 
function y?(£,t) is completely determined at time t. After that, the procedure may be 
repeated for the next time step. 
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